/*
Copyright (c) [2019年5月1日] [吴超]
[MBT_studio] is licensed under Mulan PSL v2.
You can use this software according to the terms and conditions of the Mulan PSL v2.
You may obtain a copy of Mulan PSL v2 at:
		 http://license.coscl.org.cn/MulanPSL2
THIS SOFTWARE IS PROVIDED ON AN "AS IS" BASIS, WITHOUT WARRANTIES OF ANY KIND,
EITHER EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO NON-INFRINGEMENT,
MERCHANTABILITY OR FIT FOR A PARTICULAR PURPOSE.
See the Mulan PSL v2 for more details.
*/
#include "stdafx.h"
#include "S_PhysicsFLIP.h"

#include <ppl.h>
#include <list>

#include "物体/S_物体.h"



#define DEF_VEC_UI_F(x, y, z) {uint32(x), uint32(y), uint32(z)}




S_PhysicsFLIP::S_PhysicsFLIP() {
	//m_Core = f_flip_创建流体单元()
	m_IdType = E_物理引擎::E_物理引擎_FLIP;

	m_迭代次数 = 4;

}

S_PhysicsFLIP::~S_PhysicsFLIP() {

}

void S_PhysicsFLIP::f_Update() {
	
	f_计算散度(1, m_Core->m_方向X0, m_Core->m_方向X, m_Core->m_粘度);
	f_计算散度(2, m_Core->m_方向X0, m_Core->m_方向X, m_Core->m_粘度);
	f_计算散度(3, m_Core->m_方向X0, m_Core->m_方向X, m_Core->m_粘度);

	f_投影(m_Core->m_方向X0, m_Core->m_方向Y0, m_Core->m_方向Z0, m_Core->m_方向X, m_Core->m_方向Y);

	f_流动(1, m_Core->m_方向X, m_Core->m_方向Y, m_Core->m_方向X0, m_Core->m_方向Y0, m_Core->m_方向Z0);
	f_流动(2, m_Core->m_方向X, m_Core->m_方向Y, m_Core->m_方向X0, m_Core->m_方向Y0, m_Core->m_方向Z0);
	f_流动(3, m_Core->m_方向X, m_Core->m_方向Y, m_Core->m_方向X0, m_Core->m_方向Y0, m_Core->m_方向Z0);

	f_投影(m_Core->m_方向X, m_Core->m_方向Y, m_Core->m_方向Z, m_Core->m_方向X0, m_Core->m_方向Y0);


	f_计算散度(0, m_Core->m_速度, m_Core->m_密度, m_Core->m_扩散量);
	f_流动(0, m_Core->m_密度, m_Core->m_速度, m_Core->m_方向X, m_Core->m_方向Y, m_Core->m_方向Z);
}

void S_PhysicsFLIP::f_应用物理变换() {
	
		
	
}

void S_PhysicsFLIP::f_取更新物体(std::vector<Mat44f>& mat, std::vector<void*>& userData) {
	

}

void S_PhysicsFLIP::f_取碰撞物体(std::vector<S_碰撞信息>& objs, uint32 scenekey) {
	

}


void S_PhysicsFLIP::f_取更新物体(std::vector<std::vector<Mat44f>>& mat, std::vector<std::vector<void*>>& userData) {


}

void S_PhysicsFLIP::f_添加物理体(S_物理体* body, uint32 scenekey) {
	
}

void S_PhysicsFLIP::f_移除物理体(S_物理体* body, uint32 scenekey) {
	
}

void S_PhysicsFLIP::f_添加物理体碰撞回调(S_物理体* body, uint32 scenekey) {
	
}

void S_PhysicsFLIP::f_计算散度(uint8 方向, float32* v, float32* v0, float32 散度) {
	const uvec3& 分辨率 = m_Core->m_域属性.m_单元属性.m_分辨率;

	float32 预计算散度 = m_步进时间 * 散度 * (分辨率.x - 2) * (分辨率.x - 2);
	//f_计算边界(b, x, x0, a, 1 + 6 * a, iter, N);
	f_线性方程求解(方向, v, v0, 预计算散度, 预计算散度*6+1);

}

void S_PhysicsFLIP::f_设置边界(uint8 方向, float32* v) {
	const uvec3& 分辨率 = m_Core->m_域属性.m_单元属性.m_分辨率;

	uvec3 dim = 分辨率;

	Concurrency::parallel_for<uint32>(0, 分辨率.y, [&](uint32 y) {
		//for (uint32 y = 0; y < dim.y; ++y) {
		Concurrency::parallel_for<uint32>(0, 分辨率.x, [&](uint32 x) {
		//for (uint32 x = 0; x < dim.x; ++x) {
			uint64 id = f_坐标转线性索引({ x, y, 0 }, dim);
			v[id] = 方向 == 3 ? -v[f_坐标转线性索引({ x, y, 1 }, dim)] : v[f_坐标转线性索引({ x, y, 1 }, dim)];

			id = f_坐标转线性索引({ x, y, dim.z - 1 }, dim);
			v[id] = 方向 == 3 ? -v[f_坐标转线性索引({ x, y, dim.z - 2 }, dim)] : v[f_坐标转线性索引({ x, y, dim.z - 2 }, dim)];
		//}
		});
		//}
	});
	

	Concurrency::parallel_for<uint32>(0, dim.z, [&](uint32 z) {
		//for (uint32 z = 0; z < dim.z; ++z) {
		Concurrency::parallel_for<uint32>(0, dim.x, [&](uint32 x) {
		//for (uint32 x = 0; x < dim.x; ++x) {
			uint64 id = f_坐标转线性索引({ x, 0, z }, dim);
			v[id] = 方向 == 2 ? -v[f_坐标转线性索引({ x, dim.y - 1, 1 }, dim)] : v[f_坐标转线性索引({ x, dim.y - 1, 1 }, dim)];

			id = f_坐标转线性索引({ x, dim.y - 1, z }, dim);
			v[id] = 方向 == 2 ? -v[f_坐标转线性索引({ x, dim.y - 2, z }, dim)] : v[f_坐标转线性索引({ x, dim.y - 2, z }, dim)];
		//}
		});
		//}
	});

	Concurrency::parallel_for<uint32>(0, dim.z, [&](uint32 z) {
	//for (uint32 z = 0; z < dim.z; ++z) {
		Concurrency::parallel_for<uint32>(0, dim.y, [&](uint32 y) {
		//for (uint32 y = 0; y < dim.y; ++y) {
			uint64 id = f_坐标转线性索引({ 0, y, z }, dim);
			v[id] = 方向 == 1 ? -v[f_坐标转线性索引({ dim.x - 1, y, z }, dim)] : v[f_坐标转线性索引({ dim.x - 1, y, z }, dim)];
			
			id = f_坐标转线性索引({ dim.x - 1, y, z }, dim);
			v[id] = 方向 == 1 ? -v[f_坐标转线性索引({ dim.x - 2, y, z }, dim)] : v[f_坐标转线性索引({ dim.x - 1, y, z }, dim)];
		//}
		});
	//}
	});


	//x
	v[f_坐标转线性索引({ 0, 0, 0 }, dim)] = (
		v[f_坐标转线性索引({ 1, 0, 0 }, dim)] + v[f_坐标转线性索引({ 0, 1, 0 }, dim)] + v[f_坐标转线性索引({ 0, 0, 1 }, dim)]
	) * 0.33f;

	//mx
	v[f_坐标转线性索引({ dim.x - 1, 0, 0 }, dim)] = (
		v[f_坐标转线性索引({ dim.x - 2, 0, 0 }, dim)] + v[f_坐标转线性索引({ dim.x - 1, 1, 0 }, dim)] + v[f_坐标转线性索引({ dim.x - 1, 0, 1 }, dim)]
	) * 0.33f;


	//y
	v[f_坐标转线性索引({ 0, dim.y - 1, 0 }, dim)] = (
		v[f_坐标转线性索引({ 1, dim.y - 1, 0 }, dim)] + v[f_坐标转线性索引({ 0, dim.y - 2, 0 }, dim)] + v[f_坐标转线性索引({ 0, dim.y - 1, 1 }, dim)]
	) * 0.33f;
	//z
	v[f_坐标转线性索引({ 0, 0, dim.z - 1 }, dim)] = (
		v[f_坐标转线性索引({ 1, 0, dim.z - 1 }, dim)] + v[f_坐标转线性索引({ 0, 1, dim.z - 1 }, dim)] + v[f_坐标转线性索引({ 0, 0, dim.z - 2 }, dim)]
	) * 0.33f;


	//mx my
	v[f_坐标转线性索引({ dim.x - 1, dim.y - 1, 0 }, dim)] = (
		v[f_坐标转线性索引({ dim.x - 2, dim.y - 1, 0 }, dim)] + v[f_坐标转线性索引({ dim.x - 1, dim.y - 2, 0 }, dim)] + v[f_坐标转线性索引({ dim.x - 1, dim.y - 1, 1 }, dim)]
	) * 0.33f;

	//mx y mz
	v[f_坐标转线性索引({ dim.x - 1, 0, dim.z - 1 }, dim)] = (
		v[f_坐标转线性索引({ dim.x - 2, 0, dim.z - 1 }, dim)] + v[f_坐标转线性索引({ dim.x - 1, 1, dim.z - 1 }, dim)] + v[f_坐标转线性索引({ dim.x - 1, 0, dim.z - 2 }, dim)]
	) * 0.33f;


	//x my mz
	v[f_坐标转线性索引({ 0, dim.y - 1, dim.z - 1 }, dim)] = (
		v[f_坐标转线性索引({ 1, dim.y - 1, dim.z - 1 }, dim)] + v[f_坐标转线性索引({ 0, dim.y - 2, dim.z - 1 }, dim)] + v[f_坐标转线性索引({ 0, dim.y - 1, dim.z - 2 }, dim)]
	) * 0.33f;

	//mx my mz
	v[f_坐标转线性索引({ dim.x - 1, dim.y - 1, dim.z - 1 }, dim)] = (
		v[f_坐标转线性索引({ dim.x - 2, dim.y - 1, dim.z - 1 }, dim)] + v[f_坐标转线性索引({ dim.x - 1, dim.y - 2, dim.z - 1 }, dim)] + v[f_坐标转线性索引({ dim.x - 1, dim.y - 1, dim.z - 2 }, dim)]
	) * 0.33f;

}


void S_PhysicsFLIP::f_计算边界(uint8 方向, float32* v, float32* v0, float32 预计算散度, float32 边界散度) {
	uvec3 dim = m_Core->m_域属性.m_单元属性.m_分辨率;

	//uvec3 dim = ivec_4转3(分辨率);

	float32 边界散度倒数 = 1.0f / 边界散度;


	//vec3* v = m_Core->m_方向;
	//vec3* v0 = m_Core->m_方向0;

	for (uint32 k = 0; k < m_迭代次数; ++k) {
		for (uint32 z = 1; z < dim.z - 1; ++z) {
			for (uint32 y = 1; y < dim.y - 1; ++y) {
				for (uint32 x = 1; x < dim.x - 1; ++x) {
					uint64 id = f_坐标转线性索引({ x, y, z }, dim);

					v[id] = v[id] + v0[f_坐标转线性索引({ x + 1, y, z }, dim)] * 边界散度倒数;
					v[id] = v[id] + v0[f_坐标转线性索引({ x - 1, y, z }, dim)] * 边界散度倒数;
					v[id] = v[id] + v0[f_坐标转线性索引({ x, y + 1, z }, dim)] * 边界散度倒数;
					v[id] = v[id] + v0[f_坐标转线性索引({ x, y - 1, z }, dim)] * 边界散度倒数;
					v[id] = v[id] + v0[f_坐标转线性索引({ x, y, z + 1 }, dim)] * 边界散度倒数;
					v[id] = v[id] + v0[f_坐标转线性索引({ x, y, z - 1 }, dim)] * 边界散度倒数;
				}
			}
		}

		f_设置边界(方向, v);
	}
}

void S_PhysicsFLIP::f_线性方程求解(uint8 方向, float32* v, float32* v0, float32 a, float32 c) {
	uvec3 dim = m_Core->m_域属性.m_单元属性.m_分辨率;

	//uvec3 dim = ivec_4转3(分辨率);

	float32 倒数 = 1.0 / c;

	for (uint32 k = 0; k < m_迭代次数; ++k) {
		for (uint32 z = 1; z < dim.z - 1; ++z) {
			for (uint32 y = 1; y < dim.y - 1; ++y) {
				for (uint32 x = 1; x < dim.x - 1; ++x) {
					uint64 id = f_坐标转线性索引({ x, y, z }, dim);
					v[id] = v0[id] + a * (
						v[f_坐标转线性索引({ x + 1, y, z }, dim)]
						- v[f_坐标转线性索引({ x - 1, y, z }, dim)]
						+ v[f_坐标转线性索引({ x, y + 1, z }, dim)]
						- v[f_坐标转线性索引({ x, y - 1, z }, dim)]
						+ v[f_坐标转线性索引({ x, y, z + 1 }, dim)]
						- v[f_坐标转线性索引({ x, y, z - 1 }, dim)]
						) * 倒数;
				}
			}
		}

		f_设置边界(方向, v);
	}
}

void S_PhysicsFLIP::f_投影(float32* vx, float32* vy, float32* vz, float32* p, float32* div) {
	uvec3 dim = m_Core->m_域属性.m_单元属性.m_分辨率;

	//uvec3 dim = ivec_4转3(分辨率);

	for (uint32 z = 1; z < dim.z - 1; ++z) {
		Concurrency::parallel_for<uint32>(1, dim.y - 1, [&](uint32 y) {
		//for (uint32 y = 1; y < dim.y - 1; ++y) {
			Concurrency::parallel_for<uint32>(1, dim.x - 1, [&](uint32 x) {
				//for (uint32 x = 1; x < dim.x - 1; ++x) {
				uint64 id = f_坐标转线性索引({ x, y, z }, dim);
				div[id] = -0.5f * (
					vx[f_坐标转线性索引({ x + 1, y, z }, dim)]
					- vx[f_坐标转线性索引({ x - 1, y, z }, dim)]
					+ vy[f_坐标转线性索引({ x, y + 1, z }, dim)]
					- vy[f_坐标转线性索引({ x, y - 1, z }, dim)]
					+ vz[f_坐标转线性索引({ x, y, z + 1 }, dim)]
					- vz[f_坐标转线性索引({ x, y, z - 1 }, dim)]
					) / dim.x;

				p[id] = 0;
				//}
			});
		//}
		});
	}

	f_设置边界(0, div);
	f_设置边界(0, p);
	f_线性方程求解(0, p, div, 1, 6);


	for (uint32 k = 0; k < m_迭代次数; ++k) {
		Concurrency::parallel_for<uint32>(1, dim.z - 1, [&](uint32 z) {
		//for (uint32 z = 1; z < dim.z - 1; ++z) {
			Concurrency::parallel_for<uint32>(1, dim.y - 1, [&](uint32 y) {
			//for (uint32 y = 1; y < dim.y - 1; ++y) {
				for (uint32 x = 1; x < dim.x - 1; ++x) {
					vx[f_坐标转线性索引({ x, y, z }, dim)] -= 0.5f * (p[f_坐标转线性索引({ x + 1, y, z }, dim)] - p[f_坐标转线性索引({ x - 1, y, z }, dim)]) * dim.x;
					vy[f_坐标转线性索引({ x, y, z }, dim)] -= 0.5f * (p[f_坐标转线性索引({ x, y + 1, z }, dim)] - p[f_坐标转线性索引({ x, y - 1, z }, dim)]) * dim.y;
					vz[f_坐标转线性索引({ x, y, z }, dim)] -= 0.5f * (p[f_坐标转线性索引({ x, y, z + 1 }, dim)] - p[f_坐标转线性索引({ x, y, z - 1 }, dim)]) * dim.z;
				}
			//}
			});
		//}
		});
	}

	f_设置边界(1, vx);
	f_设置边界(2, vy);
	f_设置边界(3, vz);
}


void S_PhysicsFLIP::f_流动(uint8 方向, float32* 密度, float32* 速度, float32* vx, float32* vy, float32* vz) {
	uvec3 dim = m_Core->m_域属性.m_单元属性.m_分辨率;

	//uvec3 dim = 分辨率;

	vec3 轴向步进时间 = (dim - 2) * m_步进时间;


	vec3 维度系数 = vec_转换(dim);
	vec3 XYZ减数 = { 1,1,1 };

	float32 增量密度x0, 增量密度x1, 增量密度y0, 增量密度y1, 增量密度z0, 增量密度z1;


	for (uint32 z = 1; z < dim.z - 1; ++z) {
		XYZ减数.y = 1;

		for (uint32 y = 1; y < dim.y - 1; ++y) {
			XYZ减数.x = 1;

			for (uint32 x = 1; x < dim.x - 1; ++x) {
				float32 x0, x1, y0, y1, z0, z1;

				uint64 id = f_坐标转线性索引({ x, y, z }, dim);

				vec3 tmpXYZ = vec3{ vx[id], vy[id], vz[id] } *轴向步进时间;
				vec3 xyz = XYZ减数 - tmpXYZ;

				if (xyz.x < 0.5f) xyz.x = 0.5f;
				if (xyz.x > 维度系数.x + 0.5f) xyz.x = 维度系数.x + 0.5f;
				x0 = floorf(x);
				x1 = x0 + 1.0f;

				if (xyz.y < 0.5f) xyz.y = 0.5f;
				if (xyz.y > 维度系数.y + 0.5f) xyz.y = 维度系数.y + 0.5f;
				y0 = floorf(y);
				y1 = y0 + 1.0f;

				if (xyz.z < 0.5f) xyz.z = 0.5f;
				if (xyz.z > 维度系数.z + 0.5f) xyz.z = 维度系数.z + 0.5f;
				z0 = floorf(z);
				z1 = z0 + 1.0f;



				增量密度x1 = x - x0;
				增量密度x0 = 1.0f - 增量密度x1;

				增量密度y1 = y - y0;
				增量密度y0 = 1.0f - 增量密度y1;

				增量密度z1 = z - z0;
				增量密度z0 = 1.0f - 增量密度z1;




				密度[id] = (
					增量密度x0 * (增量密度y0 * (增量密度z0 * 速度[f_坐标转线性索引(DEF_VEC_UI_F(x0, y0, z0), dim)] 
												+ 增量密度z1 * 速度[f_坐标转线性索引(DEF_VEC_UI_F(x0, y0, z1), dim)])
								+ (增量密度y1 * (增量密度z0 * 速度[f_坐标转线性索引(DEF_VEC_UI_F(x0, y1, z0), dim)]
												+ 增量密度z1 * 速度[f_坐标转线性索引(DEF_VEC_UI_F(x0, y1, z1), dim)]))
						)

					+ 增量密度x1 * (增量密度y0 *    (增量密度z0 * 速度[f_坐标转线性索引(DEF_VEC_UI_F(x1, y0, z0), dim)]  
												+ 增量密度z1 * 速度[f_坐标转线性索引(DEF_VEC_UI_F(x1, y0, z1), dim)])       
								+ (增量密度y1 * (增量密度z0 * 速度[f_坐标转线性索引(DEF_VEC_UI_F(x1, y1, z0), dim)] 
												+ 增量密度z1 * 速度[f_坐标转线性索引(DEF_VEC_UI_F(x1, y1, z1), dim)] ))      )
					);


				++XYZ减数.x;
			}

			++XYZ减数.y;
		}

		++XYZ减数.z;
	}

	f_设置边界(方向, 密度);
}















